! Compiles as 
!
! gfortran -fdefault-real-8 -fdefault-integer-8  -I/usr/local/fftw/fftw-3.3.3/include 
! -L/usr/local/fftw/fftw-3.3.3/lib -fopenmp -lfftw3_threads -lfftw3 -lm -m64  fftw_guru_test.f90 -o f_test.x
!
! remember to: export OMP_NUM_THREADS=4


program guru_test
  
  use, intrinsic :: iso_c_binding

  implicit none

  include 'fftw3.f03'

  !-- testting correctness --

  !integer, parameter          :: nx = 8              ! 3D array dimenrions   
  !integer, parameter          :: ny = 6
  !integer, parameter          :: nz = 4
  !integer                     :: nffts = 1

  !-- testing efficiency ----

  integer, parameter          :: nx =  512              ! 3D array dimenrions   
  integer, parameter          :: ny = 4096
  integer, parameter          :: nz = 4
  integer                     :: nffts = 100

  !---------------------------


  integer(kind=4), parameter  :: rank = 1
  integer(kind=4), parameter  :: hm_rank = 2
  type(FFTW_IODIM)            :: dims(1)
  type(FFTW_IODIM)            :: hm_dims(2)

  type(C_PTR)                 :: p1x, p1y, p1z
  type(C_PTR)                 :: p3xf, p3yf, p3zf
  type(C_PTR)                 :: p3xb, p3yb, p3zb


  integer, parameter          :: jx = 3              ! point at which 1D and 3D 
  integer, parameter          :: jy = 3              ! transforms to be compared
  integer, parameter          :: jz = 3

  real,    allocatable        ::  r(:,:,:)
  complex, allocatable        ::  u(:,:,:)        ! input 3D array
  complex, allocatable        ::  v(:,:,:)        ! out 3D array

  complex                     ::  ux(nx)             ! in-place 1D arrays 
  complex                     ::  uy(ny)             ! for comparison
  complex                     ::  uz(nz)

  integer                     ::  ix, iy, iz, n      ! running indices

  integer                     :: iret, nthreads, tid
  integer                     :: omp_get_thread_num,  omp_get_max_threads

 !-----------------------------------------------

  allocate (v(nx,ny,nz),  u(nx,ny,nz),  r(nx,ny,nz))

  !return

  !$OMP PARALLEL PRIVATE(tid)
  tid = omp_get_thread_num()
  nthreads = omp_get_max_threads()
  print *, 'Hello World from thread = ', tid, ' out of ', nthreads
  !$OMP END PARALLEL

  call dfftw_init_threads(iret)
  call dfftw_plan_with_nthreads(nthreads)


  !-- initiaization of data --

  call random_number(r)

  !--  1D transforms of pensils passing through jx,jy,jz point --

  call dfftw_plan_dft_1d(p1x, nx, ux, ux, FFTW_FORWARD, FFTW_ESTIMATE)
  call dfftw_plan_dft_1d(p1y, ny, uy, uy, FFTW_FORWARD, FFTW_ESTIMATE)
  call dfftw_plan_dft_1d(p1z, nz, uz, uz, FFTW_FORWARD, FFTW_ESTIMATE)

  !write(*,*) "1D-plans created"


 !--  1D transforms of 3D array in x-direction --

  dims(1)%n = nx
  dims(1)%is = 1
  dims(1)%os = 1

  hm_dims(2)%n  = ny
  hm_dims(2)%is = nx
  hm_dims(2)%os = nx

  hm_dims(1)%n  = nz
  hm_dims(1)%is = nx*ny
  hm_dims(1)%os = nx*ny

  p3xf = fftw_plan_guru_dft(rank, dims, hm_rank, hm_dims,  u, v, FFTW_FORWARD, FFTW_ESTIMATE);
  p3xb = fftw_plan_guru_dft(rank, dims, hm_rank, hm_dims,  u, v, FFTW_BACKWARD, FFTW_ESTIMATE);
 
  !write(*,*) "x-plans created"

 
 !--  1D transforms of 3D array in y-direction --

  dims(1)%n =  ny
  dims(1)%is = nx
  dims(1)%os = nx

  hm_dims(2)%n  = nx
  hm_dims(2)%is = 1
  hm_dims(2)%os = 1

  hm_dims(1)%n  = nz
  hm_dims(1)%is = nx*ny
  hm_dims(1)%os = nx*ny
 
  p3yf = fftw_plan_guru_dft(rank, dims, hm_rank, hm_dims,  u, v, FFTW_FORWARD, FFTW_ESTIMATE);
  p3yb = fftw_plan_guru_dft(rank, dims, hm_rank, hm_dims,  u, v, FFTW_BACKWARD, FFTW_ESTIMATE);

  !write(*,*) "y-plans created"

 !--  1D transforms of 3D array in z-direction --

  dims(1)%n =  nz
  dims(1)%is = nx*ny
  dims(1)%os = nx*ny

  hm_dims(2)%n  = nx
  hm_dims(2)%is = 1
  hm_dims(2)%os = 1

  hm_dims(1)%n  = ny
  hm_dims(1)%is = nx
  hm_dims(1)%os = nx

  p3zf = fftw_plan_guru_dft(rank, dims, hm_rank, hm_dims,  u, v, FFTW_FORWARD, FFTW_ESTIMATE);
  p3zb = fftw_plan_guru_dft(rank, dims, hm_rank, hm_dims,  u, v, FFTW_BACKWARD, FFTW_ESTIMATE);

  !write(*,*) "z-plans created"

 !-- Test of efficiency of transforms --

  u = r

  do n=1,nffts
    call dfftw_execute_dft(p3xf,u,u)
    call dfftw_execute_dft(p3yf,u,u)
    call dfftw_execute_dft(p3zf,u,u)
    call dfftw_execute_dft(p3xb,u,u)
    call dfftw_execute_dft(p3yb,u,u)
    call dfftw_execute_dft(p3zb,u,u)
    u = u / (nx*ny*nz)
  end do

  write(*,*) "Test completed."
  return

 !-- Test of correctness of transforms --

  u = r

  write(*,*); write(*,*) "Testing guru interface in x-direction"; write(*,*)
  ux = r(:,  jy, jz)
  call dfftw_execute_dft(p1x, ux, ux)
  call dfftw_execute_dft(p3xf, u, v)
  do ix = 1,nx
     write (*,*)  ux(ix),   v(ix,jy,jz)
  end do


  write(*,*); write(*,*) "Testing guru interface in y-direction"; write(*,*)
  uy = r(jx,  :, jz)
  call dfftw_execute_dft(p1y, uy, uy)
  call dfftw_execute_dft(p3yf, u, v)
  do iy = 1,ny
         write (*,*)  uy(iy),   v(jx,iy,jz)
  end do


  write(*,*); write(*,*) "Testing guru interface in z-direction"; write(*,*)
  uz = r(jz, jy,  :)
  call dfftw_execute_dft(p1z, uz, uz)
  call  dfftw_execute_dft(p3zf, u, v)
  do iz = 1,nz
         write (*,*)  uz(iz),   v(jx,jy,iz)
  end do

  write(*,*);  write(*,*) "Two columns of complex numbers should match.";  write(*,*)



end



